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ABSTRACT 



Context. Observations of small carbon-bearing molecules such as CCH, C4H, C-C3H2, and HCO in the Horsehead Nebula have shown 
these species to have higher abundances towards the edge of the source than towards the center. 

Aims. Given the determination of a wide range of values for f (s _1 ), the total ionization rate of hydrogen atoms, and the proposal of a 
column-dependent £(Nh), where N H is the total column of hydrogen nuclei, we desire to determine if the effects of ((N H ) in a single 
object with spatial variation can be observable. We chose the Horsehead Nebula because of its geometry and high density. 
Methods. We model the Horsehead Nebula as a near edge-on photon dominated region (PDR), using several choices for f, both 
constant and as a function of column. The column-dependent f functions are determined by a Monte Carlo model of cosmic ray 
penetration, using a steep power-law spectrum and accounting for ionization and magnetic field effects. We consider a case with 
low-metal elemental abundances as well as a sulfur-rich case. 

Results. We show that use of a column-dependent £(Nh) of 5 x 10~ 15 s _i at the surface and 7.5 x 10~ 16 s -1 at Ay = 10 on balance 
improves agreement between measured and theoretical molecular abundances, compared with constant values of (. 
Conclusions. 
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1. Introduction 

Ion-neutral reactions are the most important driving processes 
for gas-phase chemistry. Therefore it is important to understand 
the mechanisms by which chemical species in the interstellar 
medium become ionized, in order to have a more accurate pic- 
ture of the chemistry in various interstellar sources. Near the 
edge of dense clouds and throughout diffuse clouds, UV photons 
can provide a powerful ionizing force upon the medium, espe- 
cially if there is a sufficiently strong source of radiation nearby. 
These photons do not penetrate very far into dense clouds, de- 
creasing exponentially with the column density. Other ionizing 
agents, like X-rays, will penetrate farther into dense clouds, but 
deep within the object, high-energy (> 100 MeV) cosmic rays 
are the dominant ionizing force. 

The recent detection of unexpectedly large abundances of 
HJ, however, in an assortment of diffuse clouds has raised the 
old question as to whether the high ionization r ate needed is 
caused by a high flux of cosmic rays of < 1 GeV dMcCall et alJ 
120031 llndriolo et al.ll2007t) . Such low energy cosmic rays would 
not be expected to penetrate deeply into dense clouds, so that a 
column-dependent ionization rate due to cosmic rays might exist 
in denser sources. This question is best explored by examining 
the influence of cosmic ray ionization at different depths into a 
single object, and the Horsehead Nebula is an ideal candidate for 
such an investigation. 

The Horsehead Nebula, also called Barnard 33, is a dark neb- 
ula of size about 5' in the bright nebul a IC434. It is illuminate d 
by o~ Ori from a distance of about 30' dAnthonv-Twarogl| 1 9821) . 



The radiation field incident on the cloud is most commonly take n 
to be x = 60 in Draine units dDraindll97l lHabart et al.ll2005l) . 
and the geometry of the cloud is described as nearly "edge-on", 
meaning that the line between cr Ori and the Horsehead Nebula 
is nearly perpendicular to the line of sight. This makes the 
Horsehead Nebula ideal for observing column-dependent vari- 
able s in a single source. I t has an ambient magnetic field of < 6 
liG (Zaritsky et al. 1986) and a steep density gradient ranging 
from 10 2 to 10 s cm , and contains a pre-stellar core as well as 
at least one other dense region near the "throat" that will be able 
to be studied in greater detail by the Atacama Large Millimeter 
Array (ALMA) ^Ward-Thompson et al.ll2006h . 

High ab undances of small c arbon-b earing molecules were 
observed bv lTevssier et al.l d2004l) and bv lPetv et al.l (12005), with 
higher abundances of certain molecules (CCH, C-C3H2, C4H) 
observed near the edge than at the center. This led Pety et al. 
to posit that polycyclic aromatic hydrocarbons (PAH's) near the 
edge of the cloud are being destroyed by incident radiation, and 
that the products of their destruction are these small hydrocar- 
bons. A number of other molecul es have been detected, in clud- 
ing the ions HCO + and HOC + dGoicoechea et all l2009bl) , the 
carbon- bearing neutrals HC O dGerin et al. 1 120091) . I-C3H, and 
C-C3H dTeyssier et al.ll2004[). and t he sulfur-bearing species CS 
and HCS + dGoicoechea et al.l2006l) . although, except for HCO + , 
little information of their column dependence is available. 

Chemical modeling of t he Hor s ehead Nebula was dis cussed 
by IWinnewisser & Herbsj dl993l) . iTevssier et al.l (|2004) pro- 
vided the first detailed chemical PDR model of the Horsehead 
Nebula, using the Meudon PDR code dLe Bourlot et al.l IT9931 
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iLe Petit et a0 2002). A year later. iPetv etail (l2005h modeled the 
Horsehead Nebula with the same code, comparing the results 
with observations at three diff erent lines of sight, and incorpo- 
rating PAH's into the model. lHabart et alJ d2005l) determined 
a column-dependent temperature via thermal balance. None of 
these models is able to account for the high abundances of small 
hydrocarbons at the edge, or the HC3N abundance. 

Deuterium fractionation of HCO + ([DCO + ]/[HCO + ] ~ 

0.02) has been observed in the Cloud region at Ay * 10, and 

j 1 ~~i 

used to constrain its temperature to about 20 K (Petv et al. 2007). 
Neutr al atomic oxygen has also been detected dGoicoechea et alJ 
2009a), with hopes for Herschel's heightened resolution to pro- 
vide abundances of atomic oxygen for different regions in the 
cloud. 

The effect of a high er sulfur abundance was considered by 
Goicoechea et al] d2006l) . using the Meudon code (Le Peti t et all 
2006). |Petv et alJ d2007l) also used this code to better understand 
deuterium fractionation at the Horsehead e dge. The abundanc e 
of the negative ion C(,H was calculated by Mill ar et alJ d2007l) . 
although negative ions h ave not yet been observed in this region. 
lMorata & Herbstl J2008) developed a time-dependent PDR code, 
and first applied it to the Horsehead Nebula, with mixed success. 
This is the code we make use of in this paper, in tandem with 
the Meudon PDR code, which we use to determine some of the 
physical conditions. 

ICompiegne et alJ d2007h and iGoicoechea et all |2P09b) 
have performed some re cent modeling of the Horsehead re- 
gion: | Compiegne et alJ (20071) explored the dust emission. 
IGoicoechea et alJ d2009bl) self-consistently modeled the ob- 
served spatial distribution and line intensities with detailed 
depth-dependent predictions coupled with a nonlocal radiative 
transfer calculation for H 13 CO + , DC O + and HOC + . The y com- 
pared their model results with the iGerin et a D J2009V obser- 
vations of HCO + in order to constrain the electron fraction. 
Goicoechea et al. determined a very steep relative electron abun- 
dance of n e /n H ~ 10~ 4 - 10~ 8 (where n H = n(H) + 2n(H 2 )) at 
Ay ~ 0.6 - 2.0 from the cloud edge, based on a faint emission 
line attributed to HCO + near the edge. 

In this paper, we report our investigation on the effect of a 
column-dependent cosmic ray ionization rate £(Nh) on model 
results for molecular abundances and their spatial variation in 
the Horsehead Nebula. This is offered as a partial explanation 
of the high abundances of small hydrocarbons at the edge of the 
Horsehead nebula. In Section|2l we discuss the determination of 
three different £(Nh) functions, including a discussion of the role 
played by the magnetic field. In Section[3j we provide a detailed 
description of the PDR model used, and compare our calculated 
abundances with observational values using two different sets 
of elemental abundances. We also provide predicted abundances 
for observable species. In Section[4] we discuss the implications 
of these results, and a better determination of £(Nh) from single 
sources after the advent of ALMA. 



2. The Determination of £(N H ) 

The cosmic ray ionization of the interstellar medium is caused 
primarily by relativistic protons, alpha particles, and electrons. 
This ionization rate, labeled f , is typically represented as a per- 
second rate at which cosmic rays ionize atomic hydrogen. Given 
the process 

H + CR -> H + + e ~ + CR, 



Table 1 : Some values of £#2 use d m previous models 



faz (xlO " s ') 


Source 


1 nn 


Solomon & Werner ( 1 y 1 1 ) 


1 


Herbst & Klemperer (19731 


3 - 2000 


Hartauist et al. (19781 


10 - 100 


McCall et al. (2003) 


25 


Le Petit et al. (20041 


100 


Goto et al. (20081 


6-24 


Neufeld et al. (2010) 


5000 


Gupta etal. (20101 



where CR represents ionizing cosmic rays, £ is defined by the 
kinetic equation 

where the brackets signify concentration. The ionization rate of 
other species, such as H2 and He, is usually determined in chem- 
ical networks by multiplying ( by a constant. Even near the edge 
of dense clouds, the majority of hydrogen is molecular in nature, 
so it is important to note tha t, to a good approximation, (hi ~ 2£ 
dGlassgold & Langerlll974l) . 

In the last decade, results from diffuse sour ces dMcCall et al.l 

120031; ILe Petit et al1l2004t llndriolo et all2007l). including recent 
obser vations with Herschel (Gerin et al. 2010; Neufel d et al.l 
1201 01) . have most often indicated that in these environments ( 
is more than an order of magnitude higher than the generally 
accepted value of 10~ 17 s _1 . Earlier values for ( ranging from 
IP" 17 - IP" 15 s" 1 had been proposed dSpitzer & Tomaskoli 19681: 
lHartquist et al.lfl979"l:lDalgarnol2006l) . Table Q]contains a limited 
historical overview of some of the values of ( utilized in previ- 
ous models. These actually refer to molecular rather than atomic 
hydrogen. 

The observations indicating a high (, along with this wide 
range of values, led us to initiate a calculati on of column- 
depen dent functions of At the same time, iPadovani et al.l 
(2009) undertook similar calculations. They used the ionization 
and energy loss cross sections for collisi ons between cosmic 
rays a nd atomic and molecular hydrogen dCravens & Dal garno 
1978) as well as Helium to follow the flux-spectra of cosmic 
rays through a cloud and, from the flux spectra as a function of 
position, obtained the column-dependent cosmic -ray ionization 
rate for a number of initial flux-spectra. Here we report similar 
calculations but with a Monte Carlo approach in which we also 
include magnetic field effects. 



2. 1 . Initial Spectrum 

We begin by considering the form of the initial cosmic ray flux- 
spectrum j(E) (crrT 2 s _I sr~' GeV 1 per nucleon), as a function 
of energy. The spectrum has only been directly observed within 
our solar system, where the solar win d would have depleted the 
low energy cosmic rays (Parkerl ll958l) . 

Different cosmic ray spectra have been proposed based on 
assumptions about the origin of the cosmic rays. Supernova 
shocks are c urrently the favored explanation for t he origin of 
cosmic rays dBiermann et al.ll2010b lAxfordl 1 1 98 ll) . The spec- 
trum due to the supernova blast alone imposes a low-energy 
cutoff at about 100 MeV because of energy loss due to de- 
bris and strong m agnetic field effects (Havak awa et al.l fl96lb 
lip & Axfordll985l) . It is suspected that shocks in t he debris may 
re-accelerate some of the thermalized cosmic rays dip & Axfordl 
ll985Ulndrioloetal.ll2009l) . 
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Shock models favor a steep power law for low-energy cosmic 
rays, with a new cutoff at 1 MeV, below which most cosmic rays 
would again lose a significant fraction of their energy into the de- 
bris, and would either be reabsorbed into the remnant, or would 
travel too slowly to propagate throughout the galaxy. Alternate 
theories for cosmic ray acceleration exist, b ut these al so predict 
similar spectra for low-energy cosmic rays (Butt 2009). 

Comparison between measurements of the cosmic ray flux 
and theoretical cosmic ray spectra have been very useful. Basic 
statistics, "leaky-box" models, convection methods, and Monte 
Carlo methods have been applied to better constrain cosmic ray 
spectra, often by exami ning the elemental composition of the 
cosmic rays themselves. Strong et al. (2007) contains an excel- 
lent review of these methods. IWebberl dl998l) incorporated the 



newest results from Voyager into their Monte Carlo model, in 
order to determine the low energy spectrum better. 

Nevertheless, Voyager is still in a region where solar winds 
have a substantial effect. In fact, the farther the Voyager satellite 
travel s, the steeper the low ener gy spectrum becom es (IWebberl 
119981) . Indeed, as recently as IPutze et al.1 (1201 lh . statistical, 
Monte Carlo, and "leaky box" models have been unable to con- 
strain the low energy cosmic ray spectrum, due to a lack of direct 
measurement of low energy cosmic ray protons outside the solar 

infl uence. 

llndriolo et alj d2009l) list many of the proposed cos- 
mic r ay spectra. We consid er three representativ e spec- 
tra dHavakawa et al l 1 196 it ISpitzer & Tomaskol Il968t 
iNath & Biermannl 1 19941) . which are shown in Figure Q] 
These three spectra s pan the range of low energ y cosmic rays. 
The spectrum from ISpitzer & To masko (1968) is based on 
solar system measurements of the low energy cosmic ray flux, 
and contains the minimu m low-energy cosmic ray spectrum. 
INath & Biermannl d 19941) assume that the power-law for the 
cosmic ray spectrum at 1 GeV continues down until a hard 
cut-off at 1 MeV. Theirs is the highest published estimate of 
the low energy cosmic ray flux. We chose to use these three 
spectra in order to provide the full range of impact that different 
low energy cosm ic ray flux spectra h ave on the ionization rate. 
The spectrum of INath & Biermannl d 19941) increases the most 
steeply towards lower energies, that of Spit zer & Tomaskol 
1968b actually decrea ses towards lower energies, and that of 



Hayakawa et al.l (Il96ll) lies in the middle. 



2.2. Cross Sections 

We calculate the loss of energy by considering 10 4 cosmic ray 
protons with energies, E (in eV unless otherwise noted), dis- 
tributed according to the three spectra selected above. The parti- 
cles stream into a cloud of a number density n (cm -3 ). At each 
distance increment, the particles are each assigned a random 
number, which is compared with the probability of an ionizing 
or other inelastic collision over an incremental distance, deter- 
mined by cross-sections, cr (cm 2 ). For ionizing collisions by pro- 
tons we use the form of cr, from ISpitzer & Tomaskol ( 1968). The 
cross section (cm 2 ) for ionization of a hydrogen atom as a func- 
tion of E and the rest-energy of the proton (Ep) is given by 



ct.-h =7.63 x 10~ 20 1 



E + Ep 



+ 1.23 x 10~ 2U log 
-5.29x 10" 21 . 



E + Ep 



- 1 



The first term is dominant for "low" energies (500 keV < E < 50 
MeV), so for E < 50 MeV, cr hH oc l/E, down to E * 500 keV, 
when Equation (Q]) ceases to be accurate. This cross-section is 
also used below for determining the cosmic ray ionization rate of 
atomic hydrogen. For molecular hydrogen, we simply multiply 
the cross section by a factor of 2. 

Inelastic collisions are considered for atomic and molecular 
hydro gen only, and we use the cross-sections from lCravens et all 
dl97l . accounting for rotational, vibrational and electronic ex- 
citation as well as dissociation reactions . For the ionization of 
helium, the differential cross-section from Dalga rno et al.l (l999) 
is integrated to yield a total cross section: 



cr, ;He = l.5e Q A(E), 



(2) 



where A(E) (oc l/E for E < 100 MeV; oc log(£) for E > 1 
GeV) and e are parameters fit to the measurements of Shah et al. 
d!987h . 



2.3. Energy Loss 

The energy loss calculation assumes a column great enough that 
the cosmic rays will collide with gaseous atoms and molecules 
many times. Since our model is one-dimensional, we do not con- 
sider the effects of elastic collisions on the exclusion of low en- 
ergy cosmic rays from molecular clouds. 

Because there are many collisions, we are justified in uti- 
lizing the average energy lost by cosmic ray in an ionizing 
collision, W. This is equal to the ionization energy plus the 
average energy of the ejected electron. For molecular hydro- 
gen, this is determined from the differential cr oss-section by 
ICravens & Dalgarnol (fl?78l) : lDalgarno et alJdl999l) to be 



W (eV) = 7.92£ 0082 + 4.76, 



(3) 



(1) 



where E (eV) is the energy of the cosmic ray before the ionizing 
event. The energy losses from other types of inelastic collisions 
with molecular hydrogen, as well as ionizing and other inelastic 
collisions with H an d He are taken from the detailed forms in 
ICravens et "all (1 19751) . 

This energy loss is subtracted from the initial energy of 
the cosmic ray, and becomes the new energy. At each incre- 
ment, a new flux-spectrum, j(E, Nh), is calculated, and new ran- 
dom numbers are assigned to the cosmic rays. Because of the 
energy-dependence of the <x functions, lower energy cosmic rays 
have more ionizing coll isions. In the case of the spectrum of 
Nath & Biermann (1994), cosmic rays with E < 50 MeV con- 
tribute 99% to the value of ( (see Section [2~3T >. To complicate 
matters, however, there is energy loss from magnetic effects in 
addition to the loss from collisions. Magnetic energy loss is as- 
signed based on interactions with Alfven waves, as discussed 
below, using a static magnetic field of 3 yuGauss. 

2.4. Magnetic Field Effects 

Magnetic fields play an important role in the transport of cos- 
mic rays. The Lorentz force is the largest magnetic force acting 
on cosmic rays, and affects energy loss by increasing the path 
length cosmic rays travel, as they spiral along the magnetic field 
lines. This resulting increase in path length is not, however, the 
primary source of energy loss. Rather, the dominant magnetic 
field effect on cosmic rays is due to irregularities in the magnetic 
field. 

Because of the neutralization of low-energy cosmic rays, 
there will be far fewer cosmic rays at the center of the cloud 
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than at the edge. Since cosmic rays are overwhelmingly posi- 
tively charged, these losses introduce a charge imbalance in the 
cloud. Electrons are attracted to the edge, and their motion gen- 
erates magnetic field irregularities moving from the center to the 
edge of the cloud with velocity va = B(AnpY 1 ^ 2 . These irregu- 
larities, called Alfven waves, are the dominant s ource of energy 
loss, as disc ussed in Skillin g~& Strong ] (Il976h. lHartauist et al] 
d 19781 1 19791) first applied the work of[Skilling & Strong! d 19761) 
to calculate cosmic ray ionization rates, and proposed different 
values of (, de pending on the obje ct. 

Following Skilling & Strong! d!976l) . we determine the 
charge imbalance using the Monte Carlo simulation with B — 0, 
and consider it in terms of a characteristic column density, A(E) 
(cm" 2 ), determined by the simulation, at which the number of 
cosmic rays will be depleted by a factor of e. This means that Nh 
must be > A(E) for cosmic rays of energy E to be significantly 
affected by magnetic field irregularities. This function will ap- 
pear later in the analysis. 

Alfven waves are driven by the charge imbalance, and are 
damped by t he friction betw een ions and the surrounding gas, as 
discussed by iMcIvorld 1975b in terms of the collision rate T (s -1 ) 
between ions and neutrals ( Dakarno &"Dickinsonl fl968h . The 
larger F is, the less effect the waves have. The static magnetic 
field enhances the damping by absorbing smaller irregularities. 
However, the larger static magnetic field also increases va, and 
thus the frequency of collisions between the cosmic rays and the 
Alfven waves. 

This mechanism for cosmic ray energy loss by Alfven wave 
effects is important for Nh < 10 24 crrT 2 when B < 6 mG and 
iih < 10 9 crrr 3 . At a given column density, cosmic ray energy 
is substantially affected by Alfven waves for energies less than 
the energy Eo- The static magnetic field outside denser regions 
is assumed to be much smaller than the field inside these re- 
gions. Because the difference between the magnetic field inside 
and outside the cloud significantly dampens the Alfven waves 
for mi d to high energy cosmic rays, Eq cannot be greater than 50 
MeV (Cesarsk v & Volklll978h . Eo is dependent on various phys- 
ical parameters of the source in question. For typical cold and 
dense interstellar conditions, n(HI) = 1 cm 4 , n# = 10 4 cm -3 , 
T = 20 K, and B — 3 uG. Under th ese conditions, the use of 
jo(E) from lNath & Biermannl dl994 leads to E = 1 MeV at 
N H = 10 19 crrr 2 , while for Nh > 10 21 cnr 2 , E = 50 MeV. 

Integrating over energies up to this cutoff value, we can 
obtain the magnetohydrodynamic solution for jic(E,Nn), the 
"In-Cloud" cosmic ray flux-spectrum at a given Nh, to be 
dSkilling&Stronglll976h 



j lc (E<E ,N H ) = 



A(E) 



E j(E ,N H ) 
ME ) 



VA p 

Nh Je>=e 



v(E') 



■In 



n 2 mvA^loNH \J 2 - 1 
j lc (E>E ,N H ) =j(E,N H ) 



(4) 
(5) 



In this expression, j{E, Nh) is the spectrum determined using the 
Monte Carlo simulation in the absence of magnetic field effects, 
the magnetic energy density Um — B 2 /2fio (erg/crrT 3 ), Q.q is 
th e gyromagnetic frequency (s _1 ), the Compton-Getting factor 
a ( Gleeson & A xfordl 19681) is 



a — —■ 



10 E dj 
9"ME)5£' 



where jo is the initial cosmic ray flux-spectrum, y = (1 - 
v 2 /c 2 y l/2 and jo = (1 - u ( 2 /c 2 )~^ 2 where vo is the velocity of a 
cosmic ray of energy Eo- 

Given a steep initial jo(E), the approximate effect of the 
Alfven waves and Lorentz Force is to shift the cosmic ray spec- 
trum, and thus the ionization rate (see next section), from ({Nh) 
to ((5Nh), so that the ionization rate decreases more strongly 
with column. 

For cosmic ray flux-spectra that are not very steep (m < 2 
for j oc E~ m ), the shift is less extreme. Of course, for the full de- 
scription of the relationship of ( to the column density, Equation 
(0 must be calculated for E < Eq. 

2.5. The Column-Dependent Ionization Rate 

The value of ({Nh) is calculated by integrating the product of the 
flux-spectrum from eq. d3J, jic{E, Nh), and o- iu (E) from eq. dTJ, 
as a function of "depth" Nh into a cloud, with various correction 
factors: 



C(N H ) = 1.8 x 



5 r 

- x 

3 Jo 



Ano- U n{E)hc{E,NH)dE. (6) 



The factor of 5/3 (ISpitzer & Tomaskol 119681: iDalgarno et ail 
1999) takes into account the additional ionization caused by sec- 
ondary electrons, while the factor of 1.8 accounts for ionization 
due to a particles (He +2 ). These particles are the second most 
important source of ionizing cosmic rays (( a « 0.8£ p ). By com- 
parison, relativistic electrons, the third most important ionizing 
source, have little effect: ( e ~ ( p /10Q. 

Three different functions for ((Nh) have been calculated, 
based on the cosmic ray flux-spectra in FigureQ] which are cho- 
sen to be widely divergent below E « 5 00 MeV to account for 
the uncertainty in the low - energy region dHayakawa etldl il 961; 
ISpitzer & Tomaskol 1 19681: lNath & Biermannl 1 1994 . Analytical 
expressions for ((Nh), used in the models below, are: 



^//,Hayaka- 



5 x 10 4 

Nh 



&Nath = 0.002(Afc)-°- 6 + 10- 17 s- 



(7) 
(8) 



and are valid for 10 24 cm 2 > ^ 5 x 10 19 cm 2 . These an- 
alytical expressions do not seem to change significantly for 100 
cm 3 < n < 10 6 cm 3 and 5 K < T < 1000 K, beyond which the 
effects of the density and temperature on magnetic field effects 
becomes significant. 

The results are depicted in Figure [2] in terms of the vi- 
sual extinction between the cloud and the UV source (Ay). 
We determine Ay » 4.3 x \0' 22 Nh, using Q efficiencies from 
Laor &Draind dl993l) with a grain distribution (in terms of the 
"radius" of the grain, a) of n oc a~ 3 5 with r m j n = 5 nm and 
'"max = 1 A'm. With these assumptions, the analytical expressions 
for ((Ay) are: 



b^.Hayaka' 



2.2 x 10~ 17 



+ 10 -17 s -1 , 



(9) 



&Nath = 3.05 x 10- I6 (A v )-°- 6 + 10- 17 B" 1 . (10) 

These expressions are later referred to as "mid-range" and 
'"high-range" values, respectively. 

The wide range of the ionization rate demonstrates the 
importance of low-energy cosmic rays, especially at low Nh 
or Ay. Other calculations of ((Nh) have been performed, ei- 
ther for high column densities (> 10 24 cm 2 ) where low 
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E (GeV) 

Fig. 1 : Three d iff erent cosmic ray flux s pectra, taken from 
lHavakawa et al.l (1 19611) (dashed line), Spitzer & Tomaskol 
( 1968}) (dotted line), and lNath & Biermannl d 19941) (solid line). 

energy cosmic rays do not penetrate ( Umeba vashi & Na kano 
ll98lHFinocchi & Gailll 1997b, with out consideration of the mag- 
netic field dPadovani et al.1 [2009]). or in regions w here there 
are no ionizatio n losse s dPadoan & Scald 120051) . Recently, 
Pado vani & Gallil (1201 lb have incorporated the effect of mag- 
netic mirroring, whereas we have treated the effects of Alfven 
waves on cosmic ray streaming. The results in this paper sug- 
gest that Alfven waves may have a more substantial effect on 
with factor of ~10 impact on f at certain Nh for Alfven waves 
versus a factor of ~2-4 impact on £ from magnetic mirroring. 
Ultimately, a robust magnetohydrodynamics simulation of cos- 
mic ray transport would be necessary to determine what mag- 
netic field effects have the most significant impact on cosmic ray 
penetration. 

In our study below, we determine the effect of four differ- 
ent functions for £ on the chemistry in the Horsehead Nebula. 
We consider the ( (Nh) funct ions based on flux spe ctra from 
iNath & Biermannl d 19941) and lHavakawa et alj d 19611) . as well 
as constant values for ( of 1CT 15 s _1 , and 10~ 17 s _1 , the latter 
of which is effectively the ((Nh) derived from the spectrum of 
ISpitzer & Tomaskol d 1968b . 

3. Modeling the Horsehead Nebula 

We have used the PDR model of Mor ata & Herbstl d2008l) with 
the OSU 03/2008 gas-phase network^ This network is a purely 
gas-phase one that treats the PDR as a semi-infinite series of 
slabs with the radiation source impinging on one edge. It does 
not account for freeze-out or any surface chemistry, aside from 
a simple approximation for H2 formation on grains and se- 
lected ion recombination processe s. Radiative transfer and self- 
shieldi ng of H2 and CO dDraine & Bertoldil [l996t iLee et alj 
1996b) are calculated in progression, starting with the slab at 
the edge. The chemistry is solved with a time-dependent gas- 
phase kinetics model for each slab. This model, like our model 
for cosmic rays, is one-dimensional (ID). Because cosmic rays 
are thought to stream in from all sides, the effects of the geom- 
etry are mostly lost in this model. However, even with cosmic 
rays streaming in from all angles, low energy cosmic rays will 
dominate at the edge, and will be absent from the center. The 

1 http://www.physics.ohio-state.edu/~eric/research.html 
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Fig. 2: The results of the one-dimensional Monte Carlo model 
for ( described in Section [2] in terms of Ay. The solid red, 
dashed gr een, and dotted b l ue lin e s derive from the flux - 
spect ra of INath & Biermannl dl994l) . lHavakawa et all dl96ll) . 
and Spitzer & Tomasko (1968), respectively. These lines fit the 
averaged result of dozens of iterations of the Monte Carlo model. 
The results from a single Monte Carlo run using the flux- 
spectrum of lHavakawa et all (Il96ll) are included (pink dotted) 
in order to show error. 



average value of ( at a slab near the edge will be close to the 
value determined from the ID Monte Carlo model. Because the 
majority of slabs near the center will not have low-energy cos- 
mic rays, the average ( near the center also be close to the ID 
value. Thus the average value of ( at different slabs of the cloud 
will be close to the ID values we use for ( found in Figure [2] 

Following iPetv et ail d2005l) . we compare, when possible, 
observations with model results for three regions at different 
optical extinctions (Ay) from the edge of the Horsehead PDR. 
These are the IR-edge (Ay = 1.56 + 0.73), IR-Peak (Ay = 
4.55 + 1.7), and the Cloud (Ay = 11.7 + 4.1). The error bars 
in Ay are based both on the beam size of the observations and 
uncertainty in the density profile of the cloud, as discussed in 
the next section. We determine the error in fractional abundance 
by taking the ratio between the observed column den sity of the 
specie s and the error in that column density, both from Petv et al.l 
J2005I) . 

3.1. Physical Conditions and Initial Chemical Abundances 

The density profiles used are taken from lHabart et al.l (|2005). 
The temperature pr ofile is calculated from thermal balance 
dLe Petit et al.| [2006). Cosmic rays heat the interstellar medium 
through the fhermalization of secondar y ele ctrons and photons, 
dField et alJI 19691 iGlassgold & Langerl(l973h . Thermal heating 
by cosmic rays begins to dominate at Ay > 3, but the thermal 
impact of different cosmic ray ionization rates is not very signif- 
icant until ( > 10~ 16 s _1 . Since even the highest ((Nh) drops to 
« 10~ 16 s _1 at the Cloud region, the temperature difference here 
between the high ((Nh) and ( — 10~ 17 s is only about 4 K. The 
density and temperature profiles are shown in Fig. [3] 

The gas density increases with spatial d istance into the neb - 
ula as a power law with an exponent /? (Hab art et al.l 12005). 
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Fig. 3: The temperature (dashed line) and density (solid line) 
profiles as functions of visua l extinction with ^w .Nath- The den- 
sity profile is in the form of lHabart et all d2005b . our equation 
TTJ, with P = L The temperature is from thermal balance 



dLe Petit et al.ll2006l) . At A v = 10, £ « 10~ 16 s -1 , which raises 



the temperature by 



4 K at the center compared to a f of 10 



which in terms of column density can be written as: 

1,8/08+1) 

n H (N H ) 



»H,0 
«H,0 



xoriH.o 



N H < N Hfi 
N H > N H ,o, 



(ID 



where /3 > 1 is a dimensionless constant used to parameterize the 
steepness of the number density, = 2 x 10 5 cm 4 , xq = 0.02 
pc is a length scale, and Nh,o = (1 + >S) _1 1 .23 x 10 22 cirT 2 is the 
column density at a depth of xq. For our analysis, we show the 
results for [3 — 1, and discuss results for both (3-1 and (3 — 4. 
The steeper density gradient impacts the UV photon flux and the 
resulting thermal balance. There are different total densities for 
the IR-edge and IR-peak regions. The difference in UV penetra- 
tion, temperature and density at different values of Ay noticeably 
impacts the chemistry. The cosmic ray ionization, however, is 
not significantly altered by the density gradient, because for the 
ranges of density of 100 to 10 5 cm -3 , f is column-dependent, but 
not density dependent. 

Other densiti es and density profiles have been proposed. 
iPetv et alJ d2005|t used several uniform number densities and 
profiles, while IGoicoechea et al.1 d2009bl) proposed a slowly 
changing piecewise function for the density, with three sections 
instead of two, reaching 2 x 10 5 cm 4 at Ay m 5 instead of 
Ay m 1 .0, as used here. Until the number density is better de- 
termined, significant uncertainties in the extinction at a given 
angular depth will persist. 

The UV radiation field impinging on the Horsehead sur- 
face has been a topic of much discussion and uncertaint y 
dAnthonv-Twarodl 1 9821: IZhou et al.lll99l lAbergel et al1l2003l) . 
Values of^ = 30 tox = 100 in Draine units dDrainef l978) have 
been proposed. We use % = 60, because this is the most com- 
monly used value for the Horsehead PDR. The external UV field 
is important to the chemistry only for the IR-edge. For the IR- 
peak and the Cloud regions, cosmic rays are the primary ionizing 
and photochemical agent. 

Th e initial chemica l abundances used for t he Horsehead 
PDR dLee et alJ Il996ah iMorata & Herbsil l2008h are listed in 



Table [2] and represent abundances for a dark cloud prior to the 
onset of a nearby star. These abundances comprise observed val- 
ues for small (less-than-six-atom) species in TMC-1 , as well as 
calculated early-time values from Smi th et al] d2004h for atoms 
and small molecules that have not been observed, based on so- 
called "low-metal" elemental abundances. 

In addition to these initial abundances, we also investigated 
cases with much higher elemental abundance s of sulfur, based 
prima rily on the analysis of CS and HCS + by IGoicoechea et aTl 
(2006), who place the total elemental sulfu r abundance with re- 
spect to « H at 3 .5 x 10~ 6 . On the other hand. lTevssier et alJ (2004) 
used a value of [S] ~ 10~ 7 , similar to the low-metal value used 
in this part of the paper. To determine the effect of raising the 
sulfur abundance, we utilized elemental abundances for sulfur, 
relative to hydrogen, of 10~ 6 and \0 5 , starting primarily from 
the neutral atomic form. 

The abundances are calculated from time t = to steady 
state (t = 5 X 10 6 yr). Since the age of the Horsehead Nebula is 
not well -determined, values fro m 10 4 - 10 6 yr have been con- 
sidered (Morata & Herbstll2008h . We focus only on the time of 
10 5 yr, because in general the calculated results are closest to 
observational values at this time. We also use this time because 
it is a reasonable age for a molecular cloud, given its size and 
velocity gradie nt (seelPound et al.ll2003l) . Time-dependence was 
investigated by IMorata & Her bst (2008) albeit with a different 
density profile from what is used here. They found that at times 
between 10 5 yr and steady-state, the abundances of carbon chain 
species in the Cloud region become sharply lower, as is found in 
standard cold dark clouds. They also investigated times as early 
as 10 4 yr, at which time the abundance profiles are flatter. Our 
calculations for carbon chain species have reached steady state 
by 10 4 years for Ay < 5. For Ay > 5, our calculations confirm 
their findings. 

In Figures|4]to|6] we show the calculated abundances of var- 
ious molecules as continuous functions of visual extinction with 
observed values in boxes to delineate the uncertainties in both 
abundance and Ay. The calculated abundances are plotted with 
two fixed values of (: 10~ 17 s _1 and 10~ 15 s _1 , as well as with 
two column-dependent ionization rates depicted in Figure[2j the 
mid-range ((Nh) (dashed green line), and the high-range £(Nh) 
(solid red line). The fixed value of £ = 10~ 17 s _I is equivalent 
to the lowest-range £,{Nh) in Figure|2l which utilizes only high- 
energy protons. Neither of the two fixed values for f is likely 
to be physically reasonable; the low value can pertain to the in- 
ner Cloud region but is less likely to pertain to a region near the 
edge, where at least some low-energy cosmic rays exist, while 
the high value is more likely to pertain only to the edge of the 
PDR. Unless specified, the low elemental abundance of sulfur is 
utilized. 



3.2. Results: C 2 H, c - C 3 H 2 and C 4 H 

Hydrocarbons are not direct tracers of (; nevertheless, an en- 
hanced ( at the surface of the Horsehead nebula may help to 
explain the high abundances of these small hydrocarbons at the 
edge. C 2 H, c - C3H2 and C4H are formed by a complex network 
of reactions, linked at least partially to the cosmic ray ionization 
rate via several sequence of reactions based on C and C + . The 
sequence involving neutral atomic C starts with the reactions: 



H 2 + CRP -> Hj + e~ + CRP 
Hj + H 2 — > Hj + H 
C + Hi -» CH + + H 2 , 



(12) 
(13) 
(14) 
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Table 2: Initial fractional abundances with respect to 
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Notes. (1) /(X) = n{X)l{n{H) + 2n(H 2 )) 

(2) The sulfur-rich cases include [S] = 10~ 6 and 10~ 5 . 



and CH + initiates a series of chemical reactions that eventu- 
ally results in C2H, c - C3H2 and C4H via recombination with 
electrons. The C + ion is produced in three ways depending on 
physical conditions: at low extinction (Ay < 2.5), it is formed 
principally by photoionization, and can reach a fractional abun- 
dance as high as 10 4 , whereas at high extinction (Ay > 4.5) it 
is formed less efficiently by the reaction between He + and CO. 
In the middle region (2 < Ay < 5), secondary photons from 
cosmic rays form a large amount of the C + . Once produced, it 
can radiatively associate with H2 to form the CH^ ion, which 
initiates a series of rea ctions similar to those initiated by CH + 
dHerbst & Millar! 12008). Because of these alternate pathways, 
small hydrocarbons may not be as sensitive to ( very close to 
the edge or deep within the Horsehead PDR. Regardless, our ro- 
bust chemical network allows us to explore in detail the effect of 
a column-dependent £ on the Horsehead Nebula. 

The model abundances for C2H, c - C3H2, and C4H vs Ay 
can be found in Figure |4] where observed abundances with es- 
timated uncertainties are plotted as boxes for the three regions 
studied: the IR-edge, the IR-peak, and the Cloud. For C2H, our 
use of temperature and density profiles seems to account for the 
observed abundance at the IR-edge, regardless of the value of (, 
probably because C2H formation is so dependent on photon ef- 
fects at the edge. The results diverge for the IR-peak and Cloud, 
where the high-range ({Nh) and ( — 10~ 15 s~' seem to do better 
than the other two choices of f . In the IR-peak, the abundances 
obtained with the high-range £(Nh) and £ = 10~ 15 s~' come 
within a factor of ~ 5 of the observed value, and are closer still 
for the Cloud region. 

For c - C3H2, and for C4H, none of the four plots comes 
particularly close to the observed values at the center of the IR- 
edge, although the curves obtained with the high-range £(Nh) 
and f = 1CT 15 s~' graze the lower portion of the observation box 
for C4H. This discrepancy suggests that, though a high surface f 
is important, there are likely other factors that must be taken into 



account, such as PAH fragmentation dPetv et al.ll2005l) . For the 
IR-peak region, the high-range £(Nh) and £ = 1CT 15 s~' models 
lead to results that graze portions of the observational boxes for 
both species , with the others models exhibiting much too low an 
abundance. Finally, for the Cloud region, the high-range ({Nh) 
and £ = 1CT 15 s _1 models do quite well for C4H, and C-C3H2. 
while the lower ionization models show reasonable agreement 
only for the latter. 

It would appear that, on balance, the results obtained with 
the high constant g and the high-range column-dependent £ are 
closer to observation in most instances for these three hydrocar- 
bons. To further distinguish between these two sets of results, 
we focus on the abundance ratios between IR-peak and Cloud 
regions for the three carbon-chain species. The ratios are taken 
at the visual extinctions where the models agree best with the ob- 
servations, and are listed in Table [3] The reason for taking these 
ratios is that we can better compare results between a fixed and 
a column-dependent ionization rate in this manner. These ratios 
are examined only as a way to distinguish between a constant 
and a column-dependent and their use beyond this function 
is severely limited. For example, the C2H emission attributed to 
the Cloud region may be from the FUV illu minated surface ( for 
an analogous example involving HCO, see Gerin et al. 2009|). It 
is likely that the observed ratios will change and will be far bet- 
ter constrained when the Horsehead Nebula is explored at higher 
angular resolution. 

For C2H and c - C3H2, the ratios are much closer to observa- 
tion for the column-dependent £(Nh) than for £ = 1CT 15 s _1 . In 
both of these cases, the ratios from the £(Nh) model are within 
a factor of 2 of the observed ratios. For f = 10~ 15 s~', model 
ratios disagree by a factor of 5-7. In the case of C4H, the ratio 
from the constant ( agrees slightly better with observations than 
for £(Nh), although the ratios of both models are close to obser- 
vation. Also, examining the C4H abundances from Figure [4] it 
is evident that, within the Cloud region, the results from £(Nh) 
are much closer to observation than the results from ( = 1CT 15 
s~'. In summary, as well as being unphysical, the results from a 
model with a constant f = 1CT 15 s~' do not agree as closely as 
the results from a model with a column-dependent £(Nh)- 

For (3 — 4 at t — 10 5 yr, the results are not significantly 
changed for the IR-edge, and the model underestimates the small 
hydrocarbon abundances for the IR-peak region by about an or- 
der of magnitude. The reasons for this seem to be as involved as 
the hydrocarbon chemistry. The most significant factor is that the 
production of these hydrocarbons at a higher density requires a 
higher f , and for /3 = 4, the density is much higher at the IR peak 
than for /3 = 1 . Also, with the steeper density gradient, photons 
are more effective at ionizing and dissociating at the IR-Edge, 
but fall off more abruptly at higher Ay. The difference in C + for- 
mation by photons between (3-1 and /3 = 4 density profiles is a 
factor of three, and only present at Ay < 2.5. 

It should finally be mentioned that thermal balance from 
photons depends somewhat on the density, and so the temper- 
ature profiles with f3 = 1 and j3 = 4 are different. At A v = 0.001, 
the gas-phase temperature for (3 — 1 is about 300 K, where for 
f3 — 4,T « 600 K. The gas-phase temperatures for the two den- 
sity profiles converge at Ay = 1, and this undoubtedly has some 
impact on the chemistry. It should be emphasized that ({Nh) is 
similar for the steep and gradual density gradients. 

3.3. Results: HC 3 N, HCO + , HCO and the electron fraction 

Only one line of the carbon-chain species HC3 N has been de- 
tected, and this with a very large beam-size dTevssier et al.l 
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Fig. 4: Fractional abundances of C2H, c - C3H2, C4H, and HC3N as functions of Ay. The boxes represent observations with error 
bars, and the lines are the model results for ( = 1CT 17 s _I (blue dashed), f = 10~ 15 s (pink dotted), and, from Figure [2j the 
mid-range ({Nh) (green dashed) and high-range £(Nh) (red solid). 



Table 3: Abundance Ratios for Carbon-chain Species Between 
IR-peak and Cloud Regions 







IR-peak/Cloud 


Species 


Obs. 


((N H ) 10-' 5 s- 1 


C 2 H 


14 


8.6 2.7 


C-C3H2 


25 


13 3.3 


C 4 H 


5 


8 4 



Notes. The model results are obtained at 10 5 yr. Those listed un- 
der £h(Nh) are for t he high-range column-dependent f in Fig. [2] 
Observations are from lPetv et alj d2005h . 



2004). We follow Teyssier's tabulated value for Ay, and treat 
the emission as originating in the IR-peak, though there is some 
uncertainty about the origin of this emission. The four models 
all under-produce the observed abundance of HC3N by a lit- 
tle less than an order of magnitude or more, as can be seen in 
Figure [4] with the models with the high-range £(Nh) and the 
fixed f = 10~ 15 s _1 coming closest. 

Cyanoacetylene (HC3N) is not as dependent as the other 
species on cosmic ray ionization for much of the range of vi- 



sual extinction. Two reactions primarily lead to its formation: 

CN + C 2 H 2 -> HC 3 N + H, 
C 3 H 2 N + + e~ -> HC 3 N + H. 

At the edge, the first reaction is directly related to f through 
C2H2, but the second reaction involves C3H2N + , the formation 
of which is not strongly dependent on <f . In the Cloud region, the 
situation is reversed: C2H2 is less dependent on and C3H2N + 
is then closely linked with cosmic ray ionization. Because of the 
two channels for HC3N we expect less dependence on Ay except 
in the middle range: 1 < Ay < 5. The results, shown in Figure 
|4] roughly bear this out. Interestingly, both the observed and cal- 
culated abundances for HC3N are much lower than the initial 
value, which is taken from the TMC-1 abundance. The discrep- 
ancy with the Cloud value, over three orders of magnitude, is 
especially large and very different from the analogous cases for 
the hydrocarbons in Figure|4] 

Figure [5] contains the observations and model results for 
HCO + and HCO. Since HCO + is optically thick, the carbon- 
13 isotopologue was used for observations. H 13 CO + was ob- 
serve d in emission at « 40" from the PDR edge dGerin et alJ 
2009), corresponding to an Ay « 10, w hich is essentially 
the Cloud region. Following the analyses of iGerin et al.l ([2009) 
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and lGoicoechea et alJ (l2009bl) . we determined the abundance of 
HCO + from H 13 CO + by assuming 12 C/ 13 C = 60. A faint emis- 
sion feature attributed to H 13 CO + was also seen at * 10" from 
the PDR edge, corresponding to an Ay ~ 2 with our density 
profile, and so lies essentially at the IR-edge. 

In the immediate neighborhood of Ay = 2, however, none of 
the models produces enough HCO + , but the increase in abun- 
dance with increasing extinction is steep and by Ay = 3, all but 
possibly the f = 10~ 17 s _1 model produce a comparable res ult to 
what is observed at the IR-edge. iGoicoechea et al. I d2009bl) did 
much better fitting the HCO + abundances at the edge by includ- 
ing PAH's. They also modeled profiles for H 13 CO + and DCO + . 

For the Cloud value, all models are in reasonable agreement 
with observation for HCO + , coming within factors of 2-5 of the 
observed abundance. The formation of HCO + by cosmic rays is 
very direct at high extinction; in regions where UV photons can- 
not penetrate, it is almost solely the product of the destruction 
channel for protonated molecular hydrogen with carbon monox- 
ide. At the IR-edge, however, the UV driven formation by the 
reactions H2 + CO + and H2O + C + dominates. In all regions, 
HCO + is destroyed mainly by recombination. 

For neutral HCO, all model results are too low by an order 
of magnitude or more at both the IR-edge and Cloud regions, 
even with the rela tively fast reaction between CH2 and O (from 
iGerinet all 120091). Our results disagree with the model results 
from lGerin et al.l(l2009h and lGoicoechea et all (|2009b) partly be- 
cause the Meudon reaction network includes a formation mecha- 
nism absent in the OSU network, the photodissociation reaction 

H 2 CO + hv -> HCO + H, (15) 

where hv re presents an externa l UV photon. This reaction is also 
discussed in IGerinet al .1 (12009). Including this reaction enhances 
the HCO abundance by a factor of 5 in the PDR, bringing the 
HCO abundance within an order of magnitude of the observed 
value. 

The ionization fraction, f(e~), is a measure of elemental 
abundances, ionization rate, density, and chemistry, as well as a 
constraint on the coupling of the magnetic field to the matter in 
the cloud. The ionization fraction from our models, as shown in 
Figure [5] ranges from ~ 10~ 4 in the PDR to ~ 10~ 8 in the Cloud 
re gion. This range of fractio ns agrees generally with the profile 
in Goi coechea et aD ((2009b, their Figure 4). Their inferred pro- 
file for the ionization fraction would favor the mid-range ((N h) 
from the cosmic ray flux-spectrum of Hava kawaet alJ ( 1961b . 

For the steeper density profile with /3 = 4 and the high £(Nh), 
our results are somewhat different. The HCO + abundances are 
not significantly changed, and the modeled HCO abundances in- 
crease by a factor of two in the IR-Edge and IR-Peak regions (at 
10 5 yr). Significantly, our calculated abundance of HC3N comes 
into good agreement with the Cloud region observation; it is a 
factor of 3 higher than the observed abundance at t — 10 s yr. 

3.4. Tabulated Abundances 

Calculated fractional abundances (with respect to «h) obtained 
with the standard elemental abundances are listed for more than 
twenty species in Table |H including both observed and unde- 
tected molecules. The calculated results are for a time of 10 5 
yr and pertain to the center po ints of the IR-edge, IR-peak, and 
Cloud regions dPetv et al.ll2005l) . for which observational results 
are also shown, when available. Some of the tabulated abun- 
dances, HOC + especially, seem to be possible tracers for the cos- 
mic ray ionization, because their fractional abundance becomes 



Table 4: Observations and model results for fractional abun- 
dances at 10 s yr. 
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CH4 (1U ) 
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A O 
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< 0.01 
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HJNL. (1U ) 
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67 




1 -7 
1 / 
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riL^JN (1U ) 
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< 0.01 
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0.5 
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^ n n 1 
< 0.01 




4 
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(1U ) 


< J 
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\j.y 
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10 


1. Qd 
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11 


HOC + (10" 12 ) 
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6 




73 




27 


OH+ (10 13 ) 




4 




15 




5 


H 2 + (10- 13 ) 




8 




33 




12 


H 3 + (10- 10 ) 




0.1 




45 




37 


CH+ (10-") 




23 




10 




0.2 


C 2 H 4 + (10- 13 ) 




1.5 




19 




12 


CS (10- 8 ) 


US' 


0.04 


4.^ 


0.04 




0.01 


HCS + (10-") 




0.07 


4.0« 


0.08 




0.1 



Notes. The IR-edge is at A v = 1.56 ± 0.73, the IR-peak is at Ay = 
4.55 ± 1 .7 and the Cloud region isatA v = 1 1.7 ±4.1. The model results 
are for the high c olumn-dependent ( from Fig.[2] Observations are from 
iPetv et aU (2005) unl ess otherwise noted . 
(fl) HCO at A v ~ 4 bv lGerin et all d2009h. 

(b) HC 3 N observed at Ay « 4 b ylTevssier et all d2004l). 

(c) Upper limit at A v < 2, from lGoicoechea et al] d2009bh. 

<rf) HCO + from H 13 CO + at A v w 2 and A v ~ 10 bv lGerin et all {2009). 
The calculated abundance at Ay = 3 is in close agreement with the 
observed edge value. 

<e) Observed at A v ~ 2 bv lGoicoechea et al.l d2009bl). 

'•^ For the IR-peak and IR-edge , froml Teyssier et al. (2004). 
(s) HCS + observed at A v ~ 4 bv lGoicoechea et al.ld2006l) . 



more dependent on the extinction when f depends on column 
density, than when ( is a constant value. 

In this table, we consider only the model with the high- 
range ((Nh), because it is evident that, at least for carbon-chain 
species, use of this column-dependent f leads generally to better 
agreement with observations than models with lower ionization, 
and it is more physical than the constant high-ionization model. 
Also, we do not include the case of the steeper density profile in 
this table. Predictions are discussed below in Section l3~6l 

3.5. The Sulfur-Rich Case 

We considered sulfur-bearing species, both with the standard 
initial abundances, and also for a sulfur-rich environment. We 
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Fig. 5: Relative abundances of HCO + , HCO and the ionization fraction as functions of Ay. The boxes represent observations with 
error bars, and the lines are the model results for £ = 10~ 17 s _1 (green dashed), £ = 10~ 15 s _1 (pink dotted), and, from Figure[2] the 
mid-range ((Nh) ( blue dotted) and high-range £(Nh) (red solid). 



found that the higher the elemental sulfur (up to a relative 
abundance of 1CT 5 ), the closer the model matches observa- 
tions for sulfur-bearing molecules. Our results and those of 
iGoicoechea et al.1 (120061) for the chemistry and radiative trans- 
fer agree very well. 

The results for the observed sulfur-bearing species CS and 
HCS + vs Ay at 10 5 yr can be found in Figure [6] as a function 
of the sulfur elemental abundance. There are two sets of curves, 
depending upon the rate coefficient for the charge-exchange re- 
action 

S + H + — > S + + H, 

which can affect the abundances of CS and HCS + at low sulfur 
abundances. T his reaction has a listed rate coefficient of 1.3 x 
10~ 9 cm 3 s" 1 dPrasad & HuntressH "l980) but a more likely value 
of 1 x 10~ 14 cm 3 s" 1 has been calculated^ 

The agreement attained by increasing the elemental abun- 
dance, [S], to 10~ 5 comes at a cost: at 10 5 yr, all the carbon- 
bearing species in this scenario are reduced by up to a factor 
of 10 except at the IR-edge. This effect is most severe in the 
Cloud region. This depletion occurs in part because the high sul- 
fur abundance destroys hydrocarbons by reactions with S + and 
also with S at higher extinctions and because of the increased 
fractional ionization. The depletion of carbon-bearing species 
worsens agreement for all observed species except HCO + , which 
is brought to within a factor of 2 of observation in the Cloud re- 
gion. 

This problem may suggest that a more realistic gas-phase 
sulfur elemental abundance for the Horsehead Nebula should lie 
somew here around 10~ 6 , in agreement with IGoicoechea et al.l 
(2006). The abundances of observed and predicted molecules 
with [S] = 10~ 6 are in Table [5] for the same species as listed 
in Table [4] Even with this intermediate sulfur abundance, the 
calculated abundances of carbon chain species in particular are 
lowered considerably compared with the corresponding values 
in Table [4] leading to worse agreement with observation. 



2 This rate has been tabulated in The 
Controlled Fusion Atomic Data Center 

(http://www-cfadc.phy.ornl.gov/astro/ps/data/cx/hydrogen/rates/cti.dat I. 



HCS + CS 




Fig. 6: Relative abundances of HCS + and CS as a function of Ay. 
The boxes are the observations with error bars, and the lines are 
the model results with [S] = 10~ 5 (red), [S] = 10~ 6 (green) and 
[S] = 7.2 x 10~ 8 (blue), all using the high column-dependent £ 
from Figure [2] The soli d lines use a rate for S + H + — » S + + H 
of 1.3 x 10~ 9 cm 3 s _I (Prasa d & Huntressl [l980) and the dashed 
lines use a rate of 1 x 10~ 14 cm 3 s _1 



3.6. Some Predictions 

A high column-dependent f brings with it implications for 
chemistry in the Horsehead PDR. This column-dependent £ 
varies from *2x 10~ 16 s -1 at the IR-edge to « 7 x 10~ 17 s -1 in 
the Cloud region and so leads to profiles distinctive from mod- 
els with fixed ionization rates, as can be seen for carbon-chain 
species in Figures [4]to [6] 

Also, other molecules are predicted to be in amounts in 
principle observable, and these are listed among the species in 
Tables [4] and [5] Because our £(Nh) produces reasonable abun- 
dances of C4H and HC3N in selected regions with a low elemen- 
tal abundance of sulfur, we would also expect to observe, al- 
beit with some difficulty, the more complex carbon-chains C6H 
and HC5N, based on our predictions for these regions. In addi- 
tion, the molecule HCN should definitely be present in observ- 
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Table 5: Observations and model results for fractional abun- 
dances with [S] = 1(T 6 at 10 5 yr. 



Species 
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<0.01 




13 


2 
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Notes. See Table|4]for footnotes. 



able quantities, especially in inner regions, and its isomer, HNC, 
should also be observed with a ratio HCN/HNC « 1. We predict 
ammonia in observable quantities at Ay > 4, for the low-sulfur 
case. 

Given the observations of high amounts of the reactive 
molecular ions OH + and H2Q + in many molecular objects 
dGerin et al.ll2010tlGupta et al.ll2010l) . it would be useful to con- 
sider predicted abundances of these species. Our model predic- 
tions for OH + , H20 + and HsO + in the Horsehead Nebula are 
contained in Tables [4] and [5] These predictions show low abun- 
dances for the first two ions that are rather independent of which 
of the three regions we consider. The basic problem is the low 
abundance of atomic hydrogen except at the border of the PDR 
dNeufeld et al.ll2010h . Even at the IR-Edge, H 3 + is more than 
an order of magnitude higher than either OH + or H20 + , though 
none of these species should be sufficiently abundant to be de- 
tected. In the Cloud Region, where the electron density is at the 
low level of a cold dark cloud, H30 + is depleted rather slowly 
by reactions with electrons, and should achieve a high enough 
column to be detectable. 

4. Discussion 

We have modeled the Horsehead Nebula as a PDR with time- 
dependent gas-phase chemistry using a column-dependent cos- 
mic ray ionization rate £(Nh), as well as the temperature and 



density profiles of Haba rt et alJ (120051) . At a cloud age of 10 5 yr, 
the incorporation of a high £(Nh) improves agreement between 
model and observation for the small carbon-bearing molecules 
HCO + , HC 3 N, C 2 H, c - C 3 H 2 , and C 4 H compared with a more 
standard constant ionization rate. With a higher abundance of 
elemental sulfur than our standard value, the results for small 
sulfur-bearing species are improved, but at the expense of our 
calculated values for carbon-chain species. There are also pre- 
dictions of abundances and profiles for other species, some not 
yet observed in the Horsehead Nebula, which should be in prin- 
ciple observable, including HCN, HNC, NH 3 , C 6 H, HC 5 N, and 
H30 + . Some of these predictions are strongly affected, however, 
by an increase in the assumed sulfur elemental abundance. 

Our results for C-C3H2 and C4H (but not for C2H) also indi- 
cate that the fracturing of PAH's may play an important role in 
the production of these molecules towards the edge of the PDR, 
but our model does not incorpo rate the effects of PAH's. Strong 
aromatic emission, observed by Compiegn e et alJ d2007l) . poses 
some problems, however, for the hypothesis that PAH fracturing 
is the source of small hydrocarbons. These authors claim a high 
concentration of neutral PAH's in the HII region, which suggests 
that PAH's may endure the radiation at the IR-edge, instead of 
breaking apart into the observed hydrocarbons. 

The detailed form of the calculated abundance profiles in 
Figures|4]through|6]cannot be observed because observations up 
to the present lack sufficient resolution, and because the density 
profile is not well-determined. With the advent of the Atacama 
Large Millimeter Array (AL MA), the estim ated increase in an- 
gular resolution, to ~ 0.1 " dWoottenll200l . should allow us 
to observe the form of these abundance profiles, so as to bet- 
ter determine the initial flux-spectrum for cosmic rays for the 
Horsehead Nebula. 

It appears, from llndriolo et alJ d2010l) . that there is some en- 
vironmental influence on the low energy flux of cosmic rays. It 
would be of great interest to not only examine the Horsehead 
Nebula at greater angular resolution, but to also observe and 
model other PDR's such as the Orion Bar, IC-63, L1688-W, and 
portions of Sgr B2 to determine how the low energy cosmic ray 
flux varies in our Galaxy. Sgr B2 is of special interest given the 
high values of <C inferred from H3 observations in this region 
dOka et al.ll2005l) . Given the strong dependence of £ on the path 
cosmic rays travel, it is very likely that the low-energy cosmic 
ray flux will be object-dependent. 
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